function elaplace

%  numerical solution of Laplaces equation on 0 < x < 1, 0 < y < 1
%  plots numerical solution and error distribution over domain
%  error determined using the exact solution
%  matrix equation solved using cgm 

% clear all previous variables and plots
clear *
clf

% set grid parameters
n=40;
m=20;

%%%%%%%%%%%%%%%%%%  find numerical and exact solutions  %%%%%%%%%%%%%%%%%%%  

% generate grid
x=linspace(0,1,N+2);
y=linspace(0,1,M+2);
h=x(2)-x(1);
k=y(2)-y(1);
lambda=h/k;

fprintf('\n  grid used:  m = %i    N =  %i \n\n',M,N)

% generate matrix and RHS of equation
A=matrix(lambda,N,M);
b=rhs(x,y,lambda,N,M);

% use CGM to solve matrix equation
v=cgm(A,b);

% transform back to u(i,j)
u=map(x,y,v,N+2,M+2);

% calculate exact solution
sol=exact(N+2,M+2,x,y);

% calculate error
error=abs(u-sol);

fprintf('\n  Max Error in Computed Solution=  %e \n\n',max(max((error))))

%%%%%%%%%%%%%%%%%%  plot solution  %%%%%%%%%%%%%%%%%%%

% get(gcf)
%set(gcf,'Position', [1255 379 569 639]);
plotsize(569, 639)

% plot error
subplot(2,1,1)
colormap(hsv)
surf(x,y,error')               %'
view(-132,20)
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')
zlabel('Error','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);

say=['Solution of Laplaces Equation Using the CGM'];
title(say,'FontSize',14,'FontWeight','bold')
	
% plot numerical solution 
subplot(2,1,2)
surf(x,y,u')               %'
view(-132,20)
axis([0 1 0 1 -3 2]);
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')
zlabel('Solution','FontSize',14,'FontWeight','bold')
%colormap('default')
set(gca,'ztick',[-2 0 2]);
set(gca,'FontSize',14);


%%%%%%%%%%%%%%%%%%  exact solution  %%%%%%%%%%%%%%%%%%%
% exact solution:  assumes gL=gR=gB=0
function sol=exact(N,M,x,y)
sol=zeros(N,M);

%  number of Fourier modes to be used
modes=200;

for nn=1:modes
	np=nn*pi; e6=exp(6)*(-1)^nn;
	an(nn)=2*(6/5)*np*(-864*e6-120*e6*np^2+4*e6*np^4-672*np^2-5*np^4-26352)/(36+np^2)^4;
%	an(nn)=an(nn)/sinh(np);
end;

for j=1:M-1
	for i=1:N
		s=0;
		for nn=1:modes
			np=nn*pi; np1=np*(y(j)-1); np2=np*(y(j)+1);
			sinh2=(exp(np1)-exp(-np2))/(1-exp(-2*np));
			s=s+an(nn)*sinh2*sin(nn*pi*x(i));
%			s=s+an(nn)*sinh(nn*pi*y(j))*sin(nn*pi*x(i));
		end;
		sol(i,j)=s;
	end;
end;
for i=1:N
	sol(i,M)=gT(x(i));
end;

%%%%%%%%%%%%%%%%%%  boundary conditions  %%%%%%%%%%%%%%%%%%%
% top BC function (y=1)
function g=gT(x)
g=x*(1-x)*(4/5-x)*exp(6*x);

% right BC function (x=1)
function g=gR(y)
g=0;

% bottom BC function (y=0)
function g=gB(x)
g=0;

% left BC function (x=0)
function g=gL(y)
g=0;

%%%%%%%%%%%%%%%%%%  RHS of equation  %%%%%%%%%%%%%%%%%%%
% function for constructing RHS
function b=rhs(x,y,alpha,N,M)
b=zeros(N*M,1);
for i=1:N
	b(i)=alpha^2*gB(x(i+1))+b(i);
	it=(M-1)*N+i;
	b(it)=alpha^2*gT(x(i+1))+b(it);
end;
for j=1:M
	it=(j-1)*N+1;
	b(it)=gL(y(j+1))+b(it);
	it=j*N;
	b(it)=gR(y(j+1))+b(it);
end;

%%%%%%%%%%%%%%%%%%  construct A  %%%%%%%%%%%%%%%%%%%
% function for creating the matrix A
function A=matrix(alpha,N,M)

nm=N*M;

% generate the diagonal elements
D = sparse(1:nm,1:nm,2*(1+alpha^2)*ones(1,nm),nm,nm);

% generate the subdiagonal elements
SD = sparse(2:nm,1:nm-1,-ones(1,nm-1),nm,nm);
for i=N+1:N:nm-1
	SD(i,i-1)=0;
end;

% generate the sub-subdiagonal elements
SSD=sparse(N+1:nm,1:nm-N,-alpha^2*ones(1,nm-N),nm,nm);

% use symmetry of A to complete construction 
A=D+SD+SD'+SSD+SSD';

%%%%%%%%%%%%%%%%%%  convert back to u(i,j)  %%%%%%%%%%%%%%%%%%%
% function for converting back to u(i,j)
function u=map(x,y,v,N,M)
u=zeros(N,M);
for ix=1:N
	u(ix,1)=gB(x(ix));
	u(ix,M)=gT(x(ix));
end;
for iy=1:M
	u(1,iy)=gL(y(iy));
	u(N,iy)=gR(y(iy));
end;
for j=2:M-1
	for i=2:N-1
		l=(j-2)*(N-2)+i-1;
		u(i,j)=v(l);
	end;
end;

%%%%%%%%%%%%%%%%%%  CGM  %%%%%%%%%%%%%%%%%%%
% function for the CGM
function x=cgm(A,b)

% set CGM parameters
tol=0.000001;

nm=length(b);
x=zeros(nm,1);

% start iteration
r=b-A*x;
d=r;
rr=dot(r,r);
error=1;
beta=0;
counter=1;
while error > tol
	counter=counter+1;
	q=A*d;
	alpha=rr/dot(d,q);
	x=x+alpha*d;
	r0=r;
	r=r-alpha*q;
	rr0=rr;
	rr=dot(r,r);
	error=norm(alpha*d,inf);
	beta=rr/rr0;
	d=r+beta*d;
end;

fprintf('\n  Number of CGM Iterations = %i  \n  Relative Iteration Error =  %e \n\n',counter,error)

% subfunction plotsize    '
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);